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ABSTRACT  OF  OBJECTIVES  AND  ACCOMPLISHMENTS; 

It  has  been  our  objective  and  that  of  several  other  groups  to  use  Slater-type 
orbitals,  STOs,  characterized  by  exp(-R)  in  problems  of  ab  initio  quantum  chemistry 
rather  than  the  currently  used  Gaussian- type  orbitals,  GTOs,  characterized  by  exp(-R2). 
We  believe  that  STOs  will  give  more  accurate  results  because  they  can  fulfill  a  cusp 
condition  and  have  correct  asymptotic  behavior.  This  was  proven  by  us  in  the  case  of 
Hj .  But  the  real  test  would  be  for  molecules  with  two  or  more  atoms  and  electrons. 

The  general  molecular  problem  requires  at  most  four-center  molecular  integrals. 
In  just  the  last  year  we  have  devised  an  accuracy  and  fast  method  for  processing  STO 
molecular  integrals;  interior  regions  are  done  by  Gauss-Legendre  integration  after  a- 
function  expansions,  and  exterior  regions  analytically  with  the  aid  of  look-up  tables  for 
the  evaluation  of  basic  integrals.  C,  E,  and  F  matrices  that  facilitate  the  calculations  are 


stored  in  memory.  All  our  procedures  have  been  checked  by  the  computer  algebra 
program  Mathematica,  which  can  overcome  all  cancellation  errors. 

The  fact  that  several  groups  have  worked  many  years  on  the  STO  integral 
problem  without  achieving  the  goal  of  a  complete  integral  package,  reveals  the  inherent 
difficulty  of  translating  analytic  formulation  into  computer  programs.  The  most 
difficult  part,  namely,  radial  integrals,  has  now  been  solved  with  the  support  of  the 
AFOSR.  Now,  it  remains  to  use  this  strategy  to  assemble  a  complete  STO  integral 
package. 

INTRODUCTION 

The  problem  of  multicenter  molecular  integrals  over  Slater-type  orbitals  (STOs)  is 
a  long-standing  one  that  has  been  pursued  in  a  systematic  manner  since  the  time  of 
Kotani,  Lowdin,  Colson,  Bartnett,  Roothaan,  and  Ruedenberg  immediately  after  World 
War  II.  Silverstone,  Harris,  Michels,  Shavitts,  and  Steinborn  continued  work  on  the 
problem  through  the  sixties  and  seventies.  Here,  at  Florida  A&M  University,  we  caused 
renewed  interest  in  the  problem  by  holding  an  international  meeting  involving  35  of  the 
most  knowledgeable  scientists  on  Slater  orbitals  (exponential-type  orbitals)  in  the 
summer  of  1981. 

Today,  several  groups  are  making  efforts  to  bring  this  problem  to  a  close: 
Steinborn  (Germany),  Rico  (Spain),  Guseinov  (Russia),  Bouferguene  (France),  Talman 
(Canada),  and  Tai  (USA).  In  spite  of  the  massive  effort,  a  general  "black  box"  integral 
subroutine  has  not  been  developed.  I  do  believe  that  we  at  Florida  A&M  University 
have  finally  penetrated  the  mystery  of  this  "intractable"  problem.  With  the  help  of  the 
computer  algebra  program  Mathematica  we  can  increase  the  precision  of  a  number 
(using  up  to  100  or  more  digits)  and  overcome  cancellation  errors.  The  knowledge  of 
the  true  value  of  an  integral  gives  us  insight  to  program  in  FORTRAN  to  achieve 
acceptable  accuracy  at  high  speed.  Our  secret  is  that  every  orbital  is  associated  with  C, 
E,  and  F  matrices  that  are  to  be  stored  in  computer  memory.  Also,  speed  is  achieved  by 
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the  extensive  use  of  look-up  tables.  To  eliminate  time-consuming  multiplications,  a 
Gauss-Legendre  numerical  quadrature  is  done  for  interior  regions.  This  strategy  should 
be  transparent  after  reading  the  attached  paper,  "Comprehensive  Strategy  for  the 
Calculation  of  Overlap  Integrals  with  Slater-Type  Orbitals."  Now  that  all  of  the  hard 
thinking  has  been  done,  it  is  just  a  matter  of  time  before  an  integral  package  can  be 
assembled. 

REVIEW  OF  PROGRESS: 

1992.  Our  study  of  STOs  began  in  1978  (Steinborn  started  in  1968).  It  is  remarkable 
that  so  much  of  our  early  work  played  a  significant  role  in  our  final  formulation  of  this 
problem.  In  retrospect,  it  appears  that  this  long  gestation  period  was  necessary. 

During  1992  we  wrote  a  paper  on  semianalytical  methods  for  four-center 
molecular  integrals.  This  paper  is  in  the  spirit  of  our  final  method.  The  potentials  were 
calculated  analytically  and  the  second  integral  was  done  numerically  by  Simpson’s  rule. 
Now  we  do  the  interior  integration  by  Gauss-Legendre  quadrature,  and  the  exterior 
regions  is  done  entirely  analytically.  For  large  values  of  parameters,  we  introduced  the 
T  and  X  matrices  to  save  computer  time.  For  small  values  of  parameters  we  use 
expansions  in  E  and  F  matrices.  In  this  paper  only  Is  orbitals  were  used  and  jmax  was 
set  to  36.  Now,  we  have  precalculated  jmax  values  for  many  different  orbitals. 

At  the  meeting,  "Current  Trends  in  Computational  Chemistry,"  Jackson, 
Mississippi,  we  presented  an  important  paper,  "Benchmarks  for  Two-Center  Exchange 
Integrals."  In  it  and  in  an  earlier  paper  on  "The  Lowdin  a-function  and  Overlap 
Integrals"  we  worked  with  the  computer  algebra  program  Maihematica.  We  showed 
how  to  easily  generate  exact  values  for  the  C,  E,  and  F  matrices.  (This  had  previously 
been  done  in  FORTRAN,  but  with  great  difficulty).  This  permitted  the  exact 
determination  for  overlap  integrals  for  the  first  time.  The  advantage  of  Maihematica  is 
that  we  can  work  at  any  degree  of  precision  and  thereby  always  overcome  cancellation 
errors.  As  a  matter  of  interest,  it  appears  that  for  realistic  problems  only  100  digits  are 
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needed  at  most.  But  our  concern  is  with  accuracy  and  speed.  This  requires  that  we  use 
FORTRAN.  However,  it  is  always  comforting  and  helpful  in  programming  to  work 
toward  the  true  answer.  The  exchange  integral  is  given  by  an  infinite  sum  of  term, 

OO 

Twelve  terms  usually  give  accuracies  to  eight  digits.  Mathematica  gives  the 

exact  value  for  each  I,.  It  may  be  that  at  a  later  date  we  have  to  use  convergence 
accelerators  as  developed  by  Weniger  or  Bouferguene. 

We  used  a  LCAO  approach  with  only  s-orbitals  to  find  the  ground  state  of  Hj. 
As  expected,  STOs  proved  to  be  more  accurate  than  GTOs.  Again,  this  was  the  case  if  a 
basis  of  STOs  with  various  angular  momenta  were  used.  Once  more,  this  was 
confirmed  for  HeH2+- 

1993.  The  one-electron  three-center  integral  is  difficult  to  converge  if  it  is  evaluated  at 
the  apex  of  an  equilateral  triangle.  With  our  new  Power  Macintosh  we  felt  we  had  a 
good  chance  to  generate  a  closed  formula  for  this  integral.  After  a  heroic  effort  we  had 
to  abandon  this  hope.  This  cleared  the  atmosphere,  and  henceforth  we  were  free  to 
employ  numerical  methods. 

The  invited  paper,  "Developments  in  Multicenter  Molecular  Integrals,"  was 
presented  in  Spain.  Gauss-Legendre  numerical  integration  was  selected  for  the  interior 
region  and  a  Gauss-Laguerre  quadrature  for  the  exterior  region. 

1994.  Once  more  an  attempt  was  made  to  obtain  closed  formulas,  but  to  no  avail.  The 
Gauss-Laguerre  integration  proved  to  be  impractical  —  too  many  quadrature  points 
were  needed  for  orbitals  with  high  quantum  numbers. 

1995.  All  possibilities  having  been  tried,  the  final  formulation  appeared  obvious  and 
inevitable.  Now,  it  was  just  a  matter  of  doing  it  first  in  Mathematica  and  then  in 
FORTRAN.  The  spirit  of  the  method  is  presented  in  the  appended  preprint, 
"Comprehensive  Strategy  for  the  Calculation  of  Overlap  Integrals  with  Slater-Type 
Orbitals."  Also  included  is  a  Mathematica  version  of  a  program  for  a  three-center 
exchange  integral. 


5 


CONCLUSION: 

We  claim  to  have  solved  the  real  difficulty  of  performing  multicenter  molecular 
integrals  over  Slater-type  orbitals  by  developing  an  accurate  and  fast  method  for  the 
radial  integrals.  The  necessary  angular  integrals  are  analytic  and  well-known.  To 
complete  the  job  of  obtaining  a  "black  box"  integral  package  with  needed  extensive  de¬ 
bugging,  it  would  be  desirable  to  work  with  another  group  that  is  committed  to  Slater- 
type  orbitals.  After  all,  many  persons  collaborated  on  Guassian-type  integrals. 
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Comprehensive  strategy  for  the  calculation  of 
overlap  integrals  with  Slater-type  orbitals 

Herbert  W.  Jones 
Department  of  Physics 
Florida  A  &  M  University 
Talleihassee,  Florida  32307 

Abstract 

A  strategy  is  presented  for  the  calculation  of  two-center  overlap 
integrals  over  Slater-type  orbitals.  Displaced  orbitals  are  expanded  in 
spherical  harmonics  with  Lowdin  a-fimctions  as  coefficients.  The 
exponentials  in  the  a-fimctions  are  expanded  leading  to  representation  in 
terms  of  stored  E  and  F  matrices.  For  a  given  precision,  the  niimber  of 
terms  needed  for  each  orbital  for  a  specified  harmonic,  and  its  displacement 
multiplied  by  its  screening  constant,  is  pre-determined  and  stored.  A  survey 
of  this  data  is  presented.  The  one-dimensional  integration  needed  for  the 
overlap  is  done  by  Gauss-Legendre  numerical  integration  over  the  interior 
region,  and  analytically  over  the  exterior.  Complete  stability  is  achieved  and 
excellent  results  obtained.  Implications  for  all  multicenter  molecular 
integrals  are  apparent. 
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L  INTRODUCTION 


Slater-type  orbitals  (STOs),  i.e.,  exponential  type  orbitals,  would  be 
desirable  for  basis  sets  in  quantum  chemistry  and  molecular  physics  because 
they  can  satisfy  the  cusp  condition  at  nuclei  and  represent  long  range 
electronic  behavior  [1].  However,  the  notorious  computational  difficulties  of 
STOs  have  led  to  their  virtual  replacement  by  Gaussian-type  orbitals  in  basis 
sets.  Nevertheless,  significant  progress  is  being  made  in  several  quarters  on 
this  long-standing  problem  of  computing  with  STOs  [2]. 

Two-center  overlap  integrals  over  STOs  are  needed  in  semiempirical 
and  ab  initio  quantum  chemical  calculations.  There  is  a  long  history, 
starting  with  Kotani  [3],  Mulliken  [3],  and  Roothaan  [3],  of  systematic 
attempts  to  obtain  accurate  and  fast  evaluations  of  these  simplest  of 
multicenter  molecrJar  integrals.  They  were  not  entirely  successful.  Even  the 
most  recent  computer  codes  [4,  5]  must  contend  with  restrictions  on 
parameter  values. 

Overlap  integrals  are  also  used  to  introduce  new  strategies  for  general 
multicenter  molecular  integrals.  Here,  we  consider  both  uses  of  the  overlap 
problem. 

One  might  say  that  all  mathematical  methods  for  solving  multicenter 
integrals  are  correct;  it  is  the  finite  word  length  of  computers  that  is  at  the 
root  of  the  problem  of  implementation.  Cancellation  errors,  i.e.,  the 
subtraction  of  nearly  equal  numbers,  are  endemic  to  computations  with 
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STOs.  These  invidious  cancellation  or  differencing  errors  often  make  results 
.  unusable.  The  problem  of  the  speed  of  computation  must  also  be  confronted 
if  STOs  are  to  be  competitive  with  GTOs.  This  consideration  has  caused  us 
to  recommend  the  use  of  a  one-dimensional  Gauss-Legendre  numerical 
quadrature  for  the  interior  region  of  space,  i.e.,  between  the  origin  and  the 
outermost  orbital  center. 

We  have  had  considerable  success  in  using  the  Lowdin  a-function 
method  for  a  limited  class  of  STOs  [6].  We  must  now  present  new  details  so 
that  this  method  can  be  efficiently  used  with  any  STO.  The  twin  goals  of 
high  speed  and  high  accuracy  have  been  achieved  here  with  overlap  integrals 
with  implications  for  all  multicenter  integrals. 

II.  THE  a-FUNCTION  EXPANSION 

A  STO  %  =  AR^  (0,  (j))  in  its  local  coordinate  system  (R,  0,  (j)) 

has  its  origin  displaced  to  (0,  0,  a)  in  the  working  coordinate  system  (r,  0,  (})). 
Its  expansion  in  spherical  harmonics  is  [7]: 


A 

'(2L  +  1XL+M)\' 

Vz 

(-If  I 

47c(f  +  M)! 

1/2 

4n(L-M)\ 

_(2i  +  lX£-M)[_ 

(1) 

where 
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2ie+My.  S3  ' 

w,-(Cflr‘-'-‘Kry-'-‘ 


O'.  y)x 


and 


r  <  a 
r>  a 


(2) 


(3) 


The  normalization  constant  A  =  (2Q  [(27/)!]~ 2  and  wrap  =  N-\-L  +  £-M, 

and  jstop  =  N  +  i .  Suzuki  [8]  has  shown  that  -M  may  be  included  in  istop. 
It  is  necessary  when  using  the  computer  algebra  program  Mathematica  [9].) 

As  pointed  out  by  Bouferguene  and  Rinaldi  [10],  the  elements  of  the  C- 
matrix  7)  grow  explosively  with  higher  harmonics  i,  which  usually 

necessitates  using  computer  algebra  with  arbitrary  accuracy  to  avoid 
intolerable  cancellation  errors.  Expansion  of  the  exponentials  in  the  a- 
function  leads  to  the  E-matrix  for  r  <  a  and  the  F-matrix  for  r  >a  [6,  11]. 
The  elements  of  these  matrices  are  small  and  not  intimidating.  Thus, 


a 


NLM 


j  max  istop 


jstop  /max 


j=0  i=0 


J-t-l 


(4) 

(5) 
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As  a  general  rule,  computation  is  faster  if  a  double  simi  is  reduced  to  a 
single  sum.  Therefore,  for  the  sake  of  simplicity,  too,  we  introduce  single 
dimensional  matrices  Y^Q)  and  Z^Q). 


^max 

(6) 

istop 

^  ,  r>  a 

j=0 

(7) 

with 


istop 

1=0 


(8) 


and 


Z/0>E/v'’“'0'.0(W 


(9) 


/=0 


These  definitions  are  slightly  different  from  the  previous  ones  because 
we  want  to  emphasize  that  a  and  r  always  are  associated  with  Now,  the 
issue  is  how  many  terms  in  the  expansions  of  the  exponentials  are  to  be 
taken,  i.e.,  the  values  of  jmax  and  imax.  In  previous  work  with  Is  orbitals 
we  simply  used  jinax  =  imax  =  36. 
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To  cover  all  possible  values  of  r  in  the  r  <  a  case  we  determine  the 
value  of  jmax  to  produce  an  error  less  than  one  part  in  10^^  of  the  exact 
value  when  r  =  a.  The  exact  value  was  determined  by  using  Mathematica 
set  to  70  digit  accuracy  using  a  C-matrix  closed  formula,  for  the  a-function. 
Fortunately,  the  values  of  Y({j)  and  also  the  terms  that  add  up  to  Z^(y) 

decrease  monotonically  after  a  few  terms,  so  there  is  no  ambiguity  in  the 
results.  But  for  some  STOs  with  very  large  quantiun  numbers  there  is  a 
slight  loss  in  accuracy  to  one  part  in  10^^.  This  could  easily  be  overcome  by 
using  double  precision,  if  need  be. 

In  Table  I  we  take  orbitals  up  to  JV  =  5,  L  =  4,  M  =  4,  ^  =  12,  and 
(^a)  =  16.  We  note  that  the  first  non-zero  vedue  of  Y(  (y)  starts  at  j  =  i 

and  increases  in  steps  of  two;  imax  is  always  equal  to  or  less  than  jmax\ 
jmax  and  imax  increase  with  (j^a)  and  i .  The  different  values  are  only 

mildly  different  from  each  other.  Therefore,  interpolation  is  hardly 
necessary;  we  may  err  on  the  side  of  caution  without  much  of  a  change  in 
imax  or  jmax  and  produce  only  a  small  increase  in  computer  time.  For 
simplicity,  we  cotald  even  just  take  imax  =  jmax.  These  values  of  jmax  and 
imax  are  to  be  stored  as  part  of  the  data  associated  with  an  orbital  as  well  as 
its  F-matrix  and  F-matrix. 
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m.  OVERLAP  INTEGRALS 


The  overlap  integral  between  two  orbitals  and  %  is  given  by 


(10) 


The  orbital  is  placed  at  the  origin  and  %  is  placed  along  the  z-axis  at 
(0,  0,  a).  Then, 

X'  =  A'  (t>)  and 

Expanding  %  about  the  origin  in  a-functions  and  performing  the 
angular  integration,  the  orthogonality  of  spherical  harmonics  dictates  that 
only  one  term  with  i  =  L'  and  M  =  M'  survives  [12].  Thus 


S^K^dr  ga,  ^r)  (11) 

0 


with 


K  = 


(-1) 


(2Q^^^\2L  +  1)  (L  +  M)!  {V  +  M)\ 
[(2N')!  W-  (2L'  +  1)  (L-M)!  {L'-M)\ 


(12) 
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The  integral  breaks  down  into  two  parts.  The  first  part  SI  is  the 
integration  from  0  to  a.  With  an  eye  to  future  developments,  we  decide  to 
use  a  24-point  gauss-Legendre  quadrature  over  this  interior  region.  We  have 
found  that  the  a-functions  are  relatively  smooth  and  therefore  numerical 
integrations  over  them  is  very  successful.  The  a-fimction  is  accurately 
evaluated  using  the  Y(  matrix  with  one  of  the  jmax  values.  If  we  tried  an 

analytical  approach  using  the  C-matrix,  we  would  be  confronted  with  terms 
of  the  form  exp[— —  Qr] .  Such  terms,  as  many  investigators  have 

discovered,  are  a  source  of  instability.  When  one  examines  exchange 
integTEils  using  expanded  a-functions,  the  product  of  two  infinite  series  must 
be  considered.  This  problem  cem  be  avoided  by  resorting  to  a  numerical 
quadrature. 

The  second  part  of  the  integration  S2  is  from  a  to  infinity.  Here 
analytic  methods  may  safely  be  employed  as  edl  exponential  terms  are  of  the 
form  exp[-(^'-hQ/'].  Thus: 


N+l 

S2  =  K'^Z,ii)^dr  (13) 

a 


Each  integration  term  can  readily  be  tiimed  into  a  sum  by  the  formula 
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(14) 


J  dr 


=  e 


A  n!  a"-* 


The  (J)  matrix  is  determined  by  using  the  Table  I  value  for  imax 
corresponding  to  the  appropriate  orbital,  harmonic,  and  (^a)  value. 


IV.  RESULTS 

The  strategy  described  here  for  the  evaluation  of  overlap  integrals 
proved  to  be  very  efficient  and  accurate.  All  overlaps  attempted  are  accurate 
to  at  least  10  parts,  as  confirmed  by  our  computer  algebra  program. 

Table  II  shows  our  results  in  calculating  overlap  integrals  done  by 
several  authors  [13].  The  second  line  of  each  calculation  is  the  exact  result 
as  determined  by  computer  algebra  [13b].  The  various  authors  provide  from 
6  decimal  digits  to  16.  Our  results  are  impressive.  We  achieve  machine 
accuracy  of  16  decimal  digits  for  (744 1 744)  with  p  =  2,  t  =  1/2  (^'  =  3,  ^  =  l), 

without  the  cumbersome  use  of  formulas  [13d].  All  the  other  of  our  wide- 
ranging  overlaps  give  14  or  15  decimal  digits,  except  the  last  case  with  the 
highest  quantum  nvunbers  (10,  0,  0 1 10,  9,  0)  [13e]:  that  gives  11  decimal 
digits.  Since  table  I  did  not  include  all  of  the  jmax  and  imax  values  needed, 
some  were  calculated  individually.  Eventually,  all  physically  acceptable 
orbitals  will  have  jmax  and  imax  values  as  data. 
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V.  CONCLUSION 


Our  new  procedure  for  overlap  integrals  is  so  straight  forward  and 
stable  that  we  can  expect  to  approach  machine  accuracy  on  all  but  unusual 
cases.  Even  here,  limited  use  of  double  precision  can  provide  16  decimal 
digits  for  all  physical  orbitals. 

The  implications  for  all  multicenter  molecular  integrals  are  clear: 
Gauss-Legendre  quadrature  for  interior  regions  should  be  used  after  the  a- 
functions  are  expanded  and  evaluated  on  a  fixed  grid.  The  exterior  regions 
that  will  only  have  negative  exponentials  may  safely  be  handled  analytically. 

An  ah  initio  calculation  on  a  molecule  would  require  the  reading  in  of 
E-  and  i^-matrices  for  each  orbital  of  a  basis  set  (together  with  their  jmax 
and  imax  values).  Then,  Y^^j)  and  would  be  calculated  over  a 

standard  grid.  Now,  we  would  be  ready  to  calculate  all  multicenter 
molecular  integrals  needed.  This  comprehensive  approach  should  advance 
the  use  of  STOs  in  molecular  problems. 
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Orbital 

1 

jmax 

imax 

NLM 

1 

2 

22 

22 

5  4  2 

6 

30 

25 

12 

36 

30 

4 

2 

36 

32 

6 

40 

35 

12 

46 

40 

16 

2 

70 

60 

6 

70 

63 

■  12 

72 

64 

5  4  3 

1 

3 

25 

21 

6 

28 

24 

12 

36 

30 

4 

3 

35 

31 

6 

40 

34 

12 

44 

38 

16 

3 

67 

59 

6 

70 

61 

12 

72 

64 

5  4  4 

1 

4 

24 

22 

6 

28 

24 

12 

34 

30 

4 

4 

34 

32 

6 

36 

34 

12 

44 

40 

16 

4 

68 

60 

6 

68 

60 

12 

70 

66 

Table  I.  The  values  ofjmax  and  imax  needed  for  16  decimal  digit  accuracy  of 
orbital  NLM  corresponding  to  (screening  constant  times  displacement 
distance)  smd  harmonic  i . 
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threeC 


1 


(*  July  6,1995  ENERGY  OF  REGION  1,2, 3,4, 5  ;  SUM  112345  *) 

(*  three  center  exchange  *) 

(*  see  H.W- Jones,  J.Comp.Chem. 12, 1217 (1991)  *) 

nn=l;hh=:0;inm=0;z=:1.2;a=2;  jmax=36;w=2*z;  zeta=2*z; 

gl  [n_,  a_]  :  =:Which[n>=0,  a^  (n+l)  /  (n+1)  ,n<-l,a'‘  (n+1)  /  (n+l) ,  True,  Log  [a]  ]  ; 
f  1  [n__,  a_,b_]  ;  =Which[n>=0,  (-Exp [-b*a]  /b^  (n+1)  *Sum[n!  /  (n-t)  !  *  (a*b)  ^ 
(n-t) , {t, 0,n} ] +n! /b^ (n+1) ) , 
n<“l, 

( (“b) ^ (-n-1) *Sum[ (“Exp[~b*a] * (-n-t-1) ! / (-n-1) ! 

/ (-b*a) ^ (-n-t)  +  (-n-t-1) ! / (-n-1) I / (-n-t) ! ) , 

{t,  1,  -n-1}  ] )  +  (-b)  ^  (-n-1)  /  (-n-l)  !  (ExpIntegralEi  [-b*a] 
-Log  [Abs  Ib]  ]  -EulerGanmia) , 

True, ExpIntegralEi  [-b*a] -Log  [Abs  [b]  ] -EulerGainina]  ; 
f2  [n_,a_„,b_]  :  =Which  [n>=0,  Exp  [-a*b]  /h^  (n+1)  *Sum[n!  /  (n-t)  I*  (b*a)  ^  (n-t)  , 
{t,0,n}], 
n<-l, 

-  (-b)  ^  (-n-l)  *Suin[-Exp  [-b*a]  *  (-n-t-1)  !  /  (-n-1)  I  / 

(-b*a) ^ (-n-t) , {t, 1, -n-1} ] 

- (-b) ^ (-n-1) / (-n-1) ! *ExpIntegralEi [-b*a] , 

True, -ExpIntegralEi [-b*a] ] ; 
f  lgam[n_,  r_,  z_]  ;  =  1/z^  (n+1)  *Gainma  [n+1,  0,  z*r]  //N; 
flgaml [n_,w_] : =f lgam[n, l,w] ; 

(*  activate  points  and  weights  *) 
b=Sqrt[2.];  d-2 . ; 
npoint=8; 

blim=d;  alim=,001; 
bina=  (blim-alim) /2; 
apb= (alim+blim) /2 ; 

hmax=12 ; 

mstop=nn+30+hmax+l; 

nn=Table [-bma*zz [ [i] ] +apb, {i, l,npoint} ] ; 
rp=Table [bma*zz [ [i] ] +apb, {i, 1, npoint } ] ; 
vmf lgainl=Table  [0,  {m,  l,mstop}  ,  {i,  1, npoint}  ]  ; 
vpf  lgaial=Table  [0,  {m,l, mstop }  ,  { i ,  1 , npoint }  ]  ; 

Do [  vmf  Igaml  [  [iti,  i]  ]  =f  Igaml  [m,  z*nn [  [i]  ]  ]  ; 
vpf  Igaml  [  [m,  i]  ]  =f  Igaml  [m,  z*rp  [  [i]  ]  ]  , 

{m,  1, mstop}  ,  {i,  1,  npoint}  ]  ; 

alfa[r_] :=Sum[y[ [q+1] ] *r^q, {q,h,30,2}] 

ul [r_] : =r^ (nn+1) *Sum[y [ [ j+1] ] *r^ j  *f Igaml [nn+j+h+1, z*r] ,{j,h,30,2}]; 

bias=2*hmax-nn+l+l; 
stop=bias+4nn-l; 
vf2b=Table  [0,  {i,  1,  stop}  ]  ; 

Do  [vf2b  [  [i]  ]  =  f 2  [i-bias,  a, w+zeta]  ,  {i,  1,  stop}  ]  ; 

abias=-nn+2*hmax+l+l; 

astop=abias+2nn-l; 

vf 2a=Table [0, {t, l,astop} ] ; 

Do [  vf 2a [ [t] ] =f 2 [t-abias, a, zeta] , {t , 1, astop} ] ; 

vf2a=vf2a//N; 

vf  2b=vf  2b/  /N; 


aa=Sqrt  [  (2*z)  (2nn+l)  /  (2nn)  1  ] 


threeC 


2 


Timing [ exchanges 0 , o ; 

Do[ 

yb=Table [z^ j  *Exp [-z*b] *Sum[ematrix [ [h+1, i+1, j+1] ] * (z*b) ^ (i-hh-h~l) , 

{i,  0,nn+hh+h-inm}  ] ,  {  j  ,  0^  jmax}  ]  //N; 

yd=Table [z^ j*Exp [-z*d] *Sum[ematrix [ [h+1, i+1, j+1] ] * (z*d) ^ (i-hh-h-l) , 

{i, O^nn+hh+h-mm} ] , { j , 0/ jmax} ] //N; 

zbsTable [z^ { j -h-l) *Sum[ematrix I [h+1, j+1, i+1] ] * (z*b) ^i, 

{i,h, jmax}] , {j,0/nn+h}] //N; 

zdsTable [z^ ( j -h-1) *S\mi[ematrix [ [h+1/ j+l,i+l]]* (z*d) ^i, 

{i,h, jmax}] , { j,0,nn+h}] //N; 

q2sSum[yb[  [j+1]  ]  *f  Igam  [nn+l+j +h,b,  z]  ,  {  j  ,h/ jmax^  2}  ]  //N; 

C3R=Table  [0,  {m^  0, 2nn+h}  ]  ; 

Do[C3R[ [nn+j“k+l] ] =C3R[ [nn+j-k+1] ] +zb[ [j+1] ] * (nn+j ) ! / (nn+j-k) ! /w^ (k+1) , 

{j,0^nn+h}, {k, 0,nn+j}] ; 

q3A=Exp  t-w*b]  *Sujn[C3R[  [m+1]  ]  *b^m,  {m,  0, 2nn+h}  ]  //N  ; 
il=0.0; 

Do[  r=rm[[i]]; 

If [r  <=  b/Volt=r^ (nn+l) *Sum[yb[ [j+1] ] *r^ j *f Igaml [nn+l+ j+h, z*r] ,  { j ,h, jmax, 2} ] , 
volt=q2/r^  (h+1)  +q3A/r^  (h+1)  -Exp  [-w*r]  /r^  (h+1)  *Sum[C3R[  [m+1]  ]  *r^m,  {m,  0, 2nn+h}  ]  ] 
If[r  <s  d,aval=Sum[yd[ [q+1] ] *r^q, {q^h, jmax}] , 

aval=Exp  [“Z*r]  *Siim[zd[  [q+1]  ]  *r^  (q-h-1)  /  {q,  0,nn+h}  ]  ]  ; 
il=il+bma*v7w[  [i]  ]  *r^  (nn+l)  *Exp  [-z*r]  *aval*volt ; 
r=rp[ [i] ] ; 

If  [r  <=  b,  volt=r^  (nn+l)  *Sum[yb[  [  j+1]  ]  j*f  Igaml  [nn+l+ j+h,  z*r]  /  {  j  ,h,  jmax,  2}  ] , 

volt=q2/r^  (h+1)  +q3A/r''  (h+1)  -Exp  [-w*r]  /r^  (h+1)  *Sum[C3R[  [m+1]  ]  *r^m,  {m,  0, 2nn+h}  ]  ] 
If[r  <s  d/ aval=Stim[yd[  [q+1]  ]  *r^q,  {q,h,  jmax}  ]  / 

aval=Exp [-z*r] *Sum[zd[ [q+1] ] *r^ (q-h-1) , {q#  0,nn+h} ] ] ; 
il=il+bma*ww[ [i] ] *r^ (nn+l) *Exp [-z*r] *aval*volt; 

,  {i, l^npoint}] ; 


(*  analytic  outside  *) 

w=2 . *z; 

i2sq2*Sum[zd[ [j+1] ] *f2 [nn-l-2h+ j , d, w] , { j^0,nn+h}] //N; 


ZC3R=Table [0, {t, 0^  3nn+2h} ] ; 

Do  [  t=:m+ j  ; 

ZC3R[ [t+l] ]=ZC3R[ [t+1] ]+zd[ [j+1] ]*C3R[ [m+1] ] , 

{ j / 0 , nn+h} , {m, 0 , 2nn+h} ] ; 

i3Rs-Sum[ZC3R[ [t+l] ] *f2 [nn-l-2h+t , d, 4 . *z] , {t,0,3nn+2h}]//N; 
i3A=q3A*Sum[zd[ [ j+1] ] *f2 [nn-2h-l+j ,d,2.*z]/{j, 0,nn+h} ] //N; 
i3=i3R+i3A; 

il23=aa'‘4/  (2h+l)  ^2*  (il+i2+i3)  *LegendreP  [h,  1 .  /Sqrt  [2  .  ]  ]  ; 
Print [h,"  ",il,"  "fi2,"  ",i3,"  ",il23]; 

exchange=exchange + i 1 2  3 

,{hj0,12}];  Print [exchange]  ] 


N [ exchange ,16] 

0.1099365385879246 


